Using molecular simulation to understand the skin barrier.

Skin's effectiveness as a barrier to permeation of water and other chemicals rests almost entirely in the outermost layer of the epidermis, the stratum corneum (SC), which consists of layers of corneocytes surrounded by highly organized lipid lamellae. As the only continuous path through the SC, transdermal permeation necessarily involves diffusion through these lipid layers. The role of the SC as a protective barrier is supported by its exceptional lipid composition consisting of ceramides (CERs), cholesterol (CHOL), and free fatty acids (FFAs) and the complete absence of phospholipids, which are present in most biological membranes. Molecular simulation, which provides molecular level detail of lipid configurations that can be connected with barrier function, has become a popular tool for studying SC lipid systems. We review this ever-increasing body of literature with the goals of (1) enabling the experimental skin community to understand, interpret and use the information generated from the simulations, (2) providing simulation experts with a solid background in the chemistry of SC lipids including the composition, structure and organization, and barrier function, and (3) presenting a state of the art picture of the field of SC lipid simulations, highlighting the difficulties and best practices for studying these systems, to encourage the generation of robust reproducible studies in the future. This review describes molecular simulation methodology and then critically examines results derived from simulations using atomistic and then coarse-grained models.

Multiscale Simulation of Ternary Stratum Corneum Lipid Mixtures: Effects of Cholesterol Composition.

Molecular dynamics simulations of mixtures of the ceramide nonhydroxy-sphingosine (NS), cholesterol, and a free fatty acid are performed to gain molecular-level understanding of the structure of the lipids found in the stratum corneum layer of skin. A new coarse-grained force field for cholesterol was developed using the multistate iterative Boltzmann inversion (MS-IBI) method. The coarse-grained cholesterol force field is compatible with previously developed coarse-grained force fields for ceramide NS, free fatty acids, and water and validated against atomistic simulations of these lipids using the CHARMM force field. Self-assembly simulations of multilayer structures using these coarse-grained force fields are performed, revealing that a large fraction of the ceramides adopt extended conformations, which cannot occur in the single bilayer in water structures typically studied using molecular simulation. Cholesterol fluidizes the membrane by promoting packing defects, and an increase in cholesterol content is found to reduce the bilayer thickness due to an increase in interdigitation of the C24 lipid tails, consistent with experimental observations. Using a reverse-mapping procedure, a self-assembled coarse-grained multilayer system is used to construct an equivalent structure with atomistic resolution. Simulations of this atomistic structure are found to closely agree with experimentally derived neutron scattering length density profiles. Significant interlayer hydrogen bonding is observed in the inner layers of the atomistic multilayer structure that are not found in the outer layers in contact with water or in equivalent bilayer structures. This work highlights the importance of simulating multilayer structures, as compared to the more commonly studied bilayer systems, to enable more appropriate comparisons with multilayer experimental membranes. These results also provide validation of the efficacy of the MS-IBI derived coarse-grained force fields and the framework for multiscale simulation.

High-throughput screening of tribological properties of monolayer films using molecular dynamics and machine learning.

Monolayer films have shown promise as a lubricating layer to reduce friction and wear of mechanical devices with separations on the nanoscale. These films have a vast design space with many tunable properties that can affect their tribological effectiveness. For example, terminal group chemistry, film composition, and backbone chemistry can all lead to films with significantly different tribological properties. This design space, however, is very difficult to explore without a combinatorial approach and an automatable, reproducible, and extensible workflow to screen for promising candidate films. Using the Molecular Simulation Design Framework (MoSDeF), a combinatorial screening study was performed to explore 9747 unique monolayer films (116 964 total simulations) and a machine learning (ML) model using a random forest regressor, an ensemble learning technique, to explore the role of terminal group chemistry and its effect on tribological effectiveness. The most promising films were found to contain small terminal groups such as cyano and ethylene. The ML model was subsequently applied to screen terminal group candidates identified from the ChEMBL small molecule library. Approximately 193 131 unique film candidates were screened with approximately a five order of magnitude speed-up in analysis compared to simulation alone. The ML model was thus able to be used as a predictive tool to greatly speed up the initial screening of promising candidate films for future simulation studies, suggesting that computational screening in combination with ML can greatly increase the throughput in combinatorial approaches to generate in silico data and then train ML models in a controlled, self-consistent fashion.

Examining the self-assembly of patchy alkane-grafted silica nanoparticles using molecular simulation.

In this work, molecular dynamics simulations are used to examine the self-assembly of anisotropically coated "patchy" nanoparticles. Specifically, we use a coarse-grained model to examine silica nanoparticles coated with alkane chains, where the poles of the grafted nanoparticle are bare, resulting in strongly attractive patches. Through a systematic screening process, the patchy nanoparticles are found to form dispersed, string-like, and aggregated phases, dependent on the combination of alkane chain length, coating chain density, and the fractional coated surface area. Correlation analysis is used to identify the ability of various particle descriptors to predict bulk phase behavior from more computationally efficient single grafted nanoparticle simulations and demonstrates that the solvent-accessible surface area of the nanoparticle core is a key predictor of bulk phase behavior. The results of this work enhance our knowledge of the phase space of patchy nanoparticles and provide a powerful approach for future screening of these materials.

Towards Molecular Simulations that are Transparent, Reproducible, Usable By Others, and Extensible (TRUE).

Systems composed of soft matter (e.g., liquids, polymers, foams, gels, colloids, and most biological materials) are ubiquitous in science and engineering, but molecular simulations of such systems pose particular computational challenges, requiring time and/or ensemble-averaged data to be collected over long simulation trajectories for property evaluation. Performing a molecular simulation of a soft matter system involves multiple steps, which have traditionally been performed by researchers in a "bespoke" fashion, resulting in many published soft matter simulations not being reproducible based on the information provided in the publications. To address the issue of reproducibility and to provide tools for computational screening, we have been developing the open-source Molecular Simulation and Design Framework (MoSDeF) software suite. In this paper, we propose a set of principles to create Transparent, Reproducible, Usable by others, and Extensible (TRUE) molecular simulations. MoSDeF facilitates the publication and dissemination of TRUE simulations by automating many of the critical steps in molecular simulation, thus enhancing their reproducibility. We provide several examples of TRUE molecular simulations: All of the steps involved in creating, running and extracting properties from the simulations are distributed on open-source platforms (within MoSDeF and on GitHub), thus meeting the definition of TRUE simulations.

Examining Tail and Headgroup Effects on Binary and Ternary Gel-Phase Lipid Bilayer Structure.

The structural properties of two- and three-component gel-phase bilayers were studied using molecular dynamics simulations. The bilayers contain distearoylphosphatidylcholine (DSPC) phospholipids mixed with alcohols and/or fatty acids of varying tail lengths, with carbon chain lengths of 12, 16, and 24 studied. Changes in both headgroup chemistry and tail length are found to affect the balance between steric repulsion and van der Waals attraction within the bilayers, manifesting in different bilayer structural properties. Lipid components are found to be located at different depths within the bilayer depending on both chain length and headgroup chemistry. The highest bilayer ordering and lowest area per tail are found in systems with medium-length tails. While longer tails can enhance van der Waals attractions, the increased tail-length asymmetry is found to induce disorder and reduce tail packing. Bulkier headgroups further increase steric repulsion, as reflected in increased component offsets and reduced tail packing. These findings help explain how bilayer composition affects the structure of gel-phase bilayers.

MoSDeF, a Python Framework Enabling Large-Scale Computational Screening of Soft Matter: Application to Chemistry-Property Relationships in Lubricating Monolayer Films.

We demonstrate how the recently developed Python-based Molecular Simulation and Design Framework (MoSDeF) can be used to perform molecular dynamics screening of functionalized monolayer films, focusing on tribological effectiveness. MoSDeF is an open-source package that allows for the programmatic construction and parametrization of soft matter systems and enables TRUE (transferable, reproducible, usable by others, and extensible) simulations. The MoSDeF-enabled screening identifies several film chemistries that simultaneously show low coefficients of friction and adhesion. We additionally develop a Python library that utilizes the RDKit cheminformatics library and the scikit-learn machine learning library that allows for the development of predictive models for the tribology of functionalized monolayer films and use this model to extract information on terminal group characteristics that most influence tribology, based on the screening data.

Influence of Single-Stranded DNA Coatings on the Interaction between Graphene Nanoflakes and Lipid Bilayers.

Using molecular dynamics simulations, it is demonstrated that a partial coating of single-stranded DNA (ssDNA) reduces the penetration depth of a graphene nanoflake (GNF) into a phospholipid bilayer by attenuating the hydrophobic force that drives the penetration. As the GNF penetrates the bilayer, the ssDNA remains adsorbed to the GNF outside of the bilayer where it shields the graphene from the surrounding water. The penetration depth is found to be controlled by the amount of ssDNA coating the GNF, with a sparser coating resulting in a deeper penetration since the ssDNA shields less of the GNF surface. As the coating density is increased, the likelihood of the GNF entering the bilayer is reduced where it instead tends to lie flat on the bilayer surface with the sugar phosphate backbone of ssDNA interacting with the hydrophilic lipid head groups. While no bilayer disruption is observed for a partially inserted ssDNA-coated GNF, a larger, bare, partially inserted GNF is found to preferentially extract phospholipids from the bilayer, offering further evidence of lipid extraction as a main cytotoxicity mechanism of GNFs. Therefore, a coating of ssDNA may reduce the cytotoxicity of GNFs by shielding the unfavorable graphene-water interaction, thus preventing graphene penetration and lipid extraction.

Investigation of the Impact of Cross-Polymerization on the Structural and Frictional Properties of Alkylsilane Monolayers Using Molecular Simulation.

Cross-linked chemisorbed n-alkylsilane (CH3(CH2)n-1Si(OH)3) monolayers on amorphous silica surfaces have been studied and their structural properties and frictional performance were compared to those of equivalent monolayers without cross-linkages. The simulations isolated for the first time the effects of both siloxane cross-linkages and the fraction of chains chemisorbed to the surface, providing insight into a longstanding fundamental question in the literature regarding molecular-level structure. The results demonstrate that both cross-linkages and the fraction of chemisorbed chains affect monolayer structure in small but measurable ways, particularly for monolayers constructed from short chains; however, these changes do not appear to have a significant impact on frictional performance.

A Transferable, Multi-Resolution Coarse-Grained Model for Amorphous Silica Nanoparticles.

Despite the ubiquity of nanoparticles in modern materials research, computational scientists are often forced to choose between simulations featuring detailed models of only a few nanoparticles or simplified models with many nanoparticles. Herein, we present a coarse-grained model for amorphous silica nanoparticles with parameters derived via potential matching to atomistic nanoparticle data, thus enabling large-scale simulations of realistic models of silica nanoparticles. Interaction parameters are optimized to match a range of nanoparticle diameters in order to increase transferability with nanoparticle size. Analytical functions are determined such that interaction parameters can be obtained for nanoparticles with arbitrary coarse-grained fidelity. The procedure is shown to be extensible to the derivation of cross-interaction parameters between coarse-grained nanoparticles and other moieties and validated for systems of grafted nanoparticles. The optimization procedure used is available as an open-source Python package and should be readily extensible to models of non-silica nanoparticles.

Composition Dependence of Water Permeation Across Multicomponent Gel-Phase Bilayers.

The permeability of multicomponent phospholipid bilayers in the gel phase is investigated via molecular dynamics simulation. The physical role of the different molecules is probed by comparing multiple mixed-component bilayers containing distearylphosphatidylcholine (DSPC) with varying amounts of either the emollient isostearyl isostearate or long-chain alcohol (dodecanol, octadecanol, or tetracosanol) molecules. Permeability is found to depend on both the tail packing density and hydrogen bonding between lipid headgroups and water. Whereas the addition of emollient or alcohol molecules to a gel-phase DSPC bilayer can increase the tail packing density, it also disturbed the hydrogen-bonding network, which in turn can increase interfacial water dynamics. These phenomena have opposing effects on bilayer permeability, which is found to depend on the balance between enhanced tail packing and decreased hydrogen bonding.

Effect of Ceramide Tail Length on the Structure of Model Stratum Corneum Lipid Bilayers.

Lipid bilayers composed of non-hydroxy sphingosine ceramide (CER NS), cholesterol (CHOL), and free fatty acids (FFAs), which are components of the human skin barrier, are studied via molecular dynamics simulations. Since mixtures of these lipids exist in dense gel phases with little molecular mobility at physiological conditions, care must be taken to ensure that the simulations become decorrelated from the initial conditions. Thus, we propose and validate an equilibration protocol based on simulated tempering, in which the simulation takes a random walk through temperature space, allowing the system to break out of metastable configurations and hence become decorrelated from its initial configuration. After validating the equilibration protocol, which we refer to as random-walk molecular dynamics, the effects of the lipid composition and ceramide tail length on bilayer properties are studied. Systems containing pure CER NS, CER NS + CHOL, and CER NS + CHOL + FFA, with the CER NS fatty acid tail length varied within each CER NS-CHOL-FFA composition, are simulated. The bilayer thickness is found to depend on the structure of the center of the bilayer, which arises as a result of the tail-length asymmetry between the lipids studied. The hydrogen bonding between the lipid headgroups and with water is found to change with the overall lipid composition, but is mostly independent of the CER fatty acid tail length. Subtle differences in the lateral packing of the lipid tails are also found as a function of CER tail length. Overall, these results provide insight into the experimentally observed trend of altered barrier properties in skin systems where there are more CERs with shorter tails present.

Molecular dynamics simulations of stratum corneum lipid mixtures: A multiscale perspective.

The lipid matrix of the stratum corneum (SC) layer of skin is essential for human survival; it acts as a barrier to prevent rapid dehydration while keeping potentially hazardous material outside the body. While the composition of the SC lipid matrix is known, the molecular-level details of its organization are difficult to infer experimentally, hindering the discovery of structure-property relationships. To this end, molecular dynamics simulations, which give molecular-level resolution, have begun to play an increasingly important role in understanding these relationships. However, most simulation studies of SC lipids have focused on preassembled bilayer configurations, which, owing to the slow dynamics of the lipids, may influence the final structure and hence the calculated properties. Self-assembled structures would avoid this dependence on the initial configuration, however, the size and length scales involved make self-assembly impractical to study with atomistic models. Here, we report on the development of coarse-grained models of SC lipids designed to study self-assembly. Building on previous work, we present the interactions between the headgroups of ceramide and free fatty acid developed using the multistate iterative Boltzmann inversion method. Validation of the new interactions is performed with simulations of preassembled bilayers and good agreement between the atomistic and coarse-grained models is found for structural properties. The self-assembly of mixtures of ceramide and free fatty acid is investigated and both bilayer and multilayer structures are found to form. This work therefore represents a necessary step in studying SC lipid systems on multiple time and length scales.

Investigating Alkylsilane Monolayer Tribology at a Single-Asperity Contact with Molecular Dynamics Simulation.

Chemisorbed monolayer films are known to possess favorable characteristics for nanoscale lubrication of micro- and nanoelectromechanical systems (MEMS/NEMS). Prior studies have shown that the friction observed for monolayer-coated surfaces features a strong dependence on the geometry of contact. Specifically, tip-like geometries have been shown to penetrate into monolayer films, inducing defects in the monolayer chains and leading to plowing mechanisms during shear, which result in higher coefficients of friction (COF) than those observed for planar geometries. In this work, we use molecular dynamics simulations to examine the tribology of model silica single-asperity contacts under shear with monolayer-coated substrates featuring various film densities. It is observed that lower monolayer densities lead to reduced COFs, in contrast to results for planar systems where COF is found to be nearly independent of monolayer density. This is attributed to a liquid-like response to shear, whereby fewer defects are imparted in monolayer chains from the asperity, and chains are easily displaced by the tip as a result of the higher free volume. This transition in the mechanism of molecular plowing suggests that liquid-like films should provide favorable lubrication at single-asperity contacts.

Perfluoropolyethers: Development of an All-Atom Force Field for Molecular Simulations and Validation with New Experimental Vapor Pressures and Liquid Densities.

A force field for perfluoropolyethers (PFPEs) based on the general optimized potentials for liquid simulations all-atom (OPLS-AA) force field has been derived in conjunction with experiments and ab initio quantum mechanical calculations. Vapor pressures and densities of two liquid PFPEs, perfluorodiglyme (CF3-O-(CF2-CF2-O)2-CF3) and perfluorotriglyme (CF3-O-(CF2-CF2-O)3-CF3), have been measured experimentally to validate the force field and increase our understanding of the physical properties of PFPEs. Force field parameters build upon those for related molecules (e.g., ethers and perfluoroalkanes) in the OPLS-AA force field, with new parameters introduced for interactions specific to PFPEs. Molecular dynamics simulations using the new force field demonstrate excellent agreement with ab initio calculations at the RHF/6-31G* level for gas-phase torsional energies (<0.5 kcal mol-1 error) and molecular structures for several PFPEs, and also accurately reproduce experimentally determined densities (<0.02 g cm-3 error) and enthalpies of vaporization derived from experimental vapor pressures (<0.3 kcal mol-1). Additional comparisons between experiment and simulation show that polyethers demonstrate a significant decrease in enthalpy of vaporization upon fluorination unlike related molecules (e.g., alkanes and alcohols). Simulation suggests this phenomenon is a result of reduced cohesion in liquid PFPEs due to a reduction in localized associations between backbone oxygen atoms and neighboring molecules.

Structural Properties of Phospholipid-based Bilayers with Long-Chain Alcohol Molecules in the Gel Phase.

The structural properties of two-component gel-phase bilayers of distearylphosphatidylcholine (DSPC) and alcohol molecules with different compositions and chain lengths (12-24 carbons long) are studied via molecular dynamics simulations. Several bilayer properties, including area per lipid, tilt angle, chain interdigitation, and headgroup offset, are studied for each system and compared, revealing important structural implications depending upon headgroup size and chain length. While tail tilt is the primary mechanism for single-component bilayers to balance tail attraction and headgroup repulsion, our results demonstrate that the lipid mixtures studied adjust this balance via an offset between the depths of the different molecular species in the bilayer; this behavior is found to depend both on composition and on the length of alcohol molecules relative to the length of DSPC tails. It is shown that the structural properties of bilayers with asymmetric tail lengths depend strongly on the bilayer composition, while the composition has less influence on mixed-component bilayers with nearly symmetric tail lengths. These findings are explained on the basis of the interdigitation between bilayer leaflets and how interdigitation is related to other structural properties.

Investigating the Structure of Multicomponent Gel-Phase Lipid Bilayers.

Single- and multicomponent lipid bilayers of 1,2-dipalmitoyl-sn-glycero-3-phosphatidylcholine (DPPC), 1,2-distearoyl-sn-glycero-3-phosphatidylcholine (DSPC), isostearyl isostearate, and heptadecanoyl heptadecanoate in the gel phase are studied via molecular dynamics simulations. It is shown that the structural properties of multicomponent bilayers can deviate strongly from the structures of their single-component counterparts. Specifically, the lipid mixtures are shown to adopt a compact packing by offsetting the positioning depths at which different lipid species are located in the bilayer. This packing mechanism affects the area per lipid, the bilayer height, and the chain tilt angles and has important consequences for other bilayer properties, such as interfacial hydrogen bonding and bilayer permeability. In particular, the simulations suggest that bilayers containing isostearyl isostearate or heptadecanoyl heptadecanoate are less permeable than pure 1,2-dipalmitoyl-sn-glycero-3-phosphatidylcholine or DSPC bilayers. Furthermore, hydrogen-bond analysis shows that the residence times of lipid-water hydrogen bonds depend strongly on the bilayer composition, with longer residence times for bilayers that have a higher DSPC content. The findings illustrate and explain the fundamental differences between the properties of single- and multicomponent bilayers.

A Coarse-Grained Model of Stratum Corneum Lipids: Free Fatty Acids and Ceramide NS.

Ceramide (CER)-based biological membranes are used both experimentally and in simulations as simplified model systems of the skin barrier. Molecular dynamics studies have generally focused on simulating preassembled structures using atomistically detailed models of CERs, which limit the system sizes and time scales that can practically be probed, rendering them ineffective for studying particular phenomena, including self-assembly into bilayer and lamellar superstructures. Here, we report on the development of a coarse-grained (CG) model for CER NS, the most abundant CER in human stratum corneum. Multistate iterative Boltzmann inversion is used to derive the intermolecular pair potentials, resulting in a force field that is applicable over a range of state points and suitable for studying ceramide self-assembly. The chosen CG mapping, which includes explicit interaction sites for hydroxyl groups, captures the directional nature of hydrogen bonding and allows for accurate predictions of several key structural properties of CER NS bilayers. Simulated wetting experiments allow the hydrophobicity of CG beads to be accurately tuned to match atomistic wetting behavior, which affects the whole system, since inaccurate hydrophobic character is found to unphysically alter the lipid packing in hydrated lamellar states. We find that CER NS can self-assemble into multilamellar structures, enabling the study of lipid systems more representative of the multilamellar lipid structures present in the skin barrier. The coarse-grained force field derived herein represents an important step in using molecular dynamics to study the human skin barrier, which gives a resolution not available through experiment alone.

Influence of Surface Morphology on the Shear-Induced Wear of Alkylsilane Monolayers: Molecular Dynamics Study.

Chemisorbed alkylsilane monolayer coatings have been shown to possess favorable lubrication properties; however, film degradation prevents the widespread use of these materials as lubricants in micro- and nanoelectromechanical systems (MEMS/NEMS). In this work, molecular dynamics (MD) simulations are used to provide insight into the conditions that promote the degradation and wear of these materials. This is achieved through removal of interfacial chain-substrate bonds during shear and the examination of the mobility of the resulting free, unbound chains. Specific focus is given to the effects of surface morphology, which has been shown previously to strongly influence frictional forces in monolayer systems. In-plane order of chain attachments is shown to lead to pressure-induced orientational ordering of monolayers, promoting film stability. This behavior is lost as nonideality is introduced into the substrate and chain patterning on the surface becomes disordered. The presence of surface roughness is found to reduce film stability, with localization of wear observed for chain attachment sites nearest the interface of contact. The influence of substrate nonideality on monolayer degradation is shown to diminish as chain length is increased.

Development of a coarse-grained water forcefield via multistate iterative Boltzmann inversion.

A coarse-grained water model is developed using multistate iterative Boltzmann inversion. Following previous work, the k-means algorithm is used to dynamically map multiple water molecules to a single coarse-grained bead, allowing the use of structure-based coarse-graining methods. The model is derived to match the bulk and interfacial properties of liquid water and improves upon previous work that used single state iterative Boltzmann inversion. The model accurately reproduces the density and structural correlations of water at 305 K and 1.0 atm, stability of a liquid droplet at 305 K, and shows little tendency to crystallize at physiological conditions. This work also illustrates several advantages of using multistate iterative Boltzmann inversion for deriving generally applicable coarse-grained forcefields.